####################################################################################################################################
# ZCTA-Level MSA Maps for COVID Paper
####################################################################################################################################

## Set globals

# In the Stata do file that called this R code, we passed two arguments:
# first, a dummy for whether the code is executing in the build
# second, the root directory location
args <- commandArgs()
build_execution <- (args[2] == "1")
root_dir <- args[3]

code_dir <- paste(root_dir, "/data/derived/R packages", sep = "")
shapefiles_dir <- paste(root_dir, "/data/derived/Shapefiles", sep = "")
data_folder <- paste(root_dir, "/data/derived/Employment", sep = "")
output_folder <- paste(root_dir, "/results/Employment", sep = "")

#***********************************************************************************************************************************
# Set-up: Load packages, set filepaths, load functions
#***********************************************************************************************************************************

# define list of packages required for this script
packages <- c("leaflet",
              "dplyr",
              "tigris",
              "leaflet.extras",
              "stringr",
              "mapview",
              "webshot",
              "scales",
              "rgeos", 
              "here")

# Load packages from local directory if we're using the runner; else, download from Internet
# i.e. we've downloaded these packages beforehand and now store them in the repo

if (build_execution == TRUE) {
  .libPaths(code_dir)
}

for (package in packages) {
  if (build_execution == TRUE) {
    library(package, character.only = TRUE, lib.loc = code_dir)
  } else {
    if (!require(package, character.only = TRUE)) {
      install.packages(package)
    }
    library(package, character.only = TRUE)
  }
}

# Install PhantomJS
webshot::install_phantomjs(force = TRUE)

# If running in build, edit webshot:::phantom_run to disable the processx::supervise dependency
if (build_execution == TRUE) {
  
  phantom_run_edited <- function(args, wait = TRUE, quiet = FALSE) {
    phantom_bin <- webshot:::find_phantom(quiet = quiet)
    
    # Handle missing phantomjs
    if (is.null(phantom_bin)) return(NULL)
    
    # Make sure args is a char vector
    args <- as.character(args)
    
    p <- callr::process$new(
      phantom_bin,
      args = args,
      stdout = "|", stderr = "|",
      supervise = FALSE
    )
    
    if (isTRUE(wait)) {
      on.exit({
        p$kill()
      })
      cat_n <- function(txt) {
        if (length(txt) > 0) {
          cat(txt, sep = "\n")
        }
      }
      while(p$is_alive()) {
        p$wait(200) # wait until min(c(time_ms, process ends))
        cat_n(p$read_error_lines())
        cat_n(p$read_output_lines())
      }
    }
    p$get_exit_status()
  }
  
  assignInNamespace("phantom_run", phantom_run_edited, ns = "webshot", pos = "package:webshot")
}

# define the dataset 
input_data <- read.csv(file.path(data_folder, "Map - Change in Low-Income Employment, ZIP.csv"),
                       stringsAsFactors = FALSE)

# Cache ZCTA spatial data so we don't have to download it every time 
tigris_cache_dir(shapefiles_dir)
options(tigris_use_cache = TRUE)

#***********************************************************************************************************************************
# Outline/structure of the rest of the script
#***********************************************************************************************************************************

## we are going to loop through cities in order to produce the final set of maps.

#***********************************************************************************************************************************
# Setting general mapping options, consistent across all maps[***]
#***********************************************************************************************************************************

# list of cities to map
cities <- c("New York City", "San Francisco", "Chicago")

#***********************************************************************************************************************************
# Looping through datasets and setting your specific options[***]
#***********************************************************************************************************************************

#-----------------------------------------------------------------------------------------------------------------------------------
# MAPPING VARIABLES & OPTIONS
#-----------------------------------------------------------------------------------------------------------------------------------

### Output extension 
file_ext <- ".png"

### MAPPING VARIABLES (n-element vector of strings for the variable names in your dataset)
map_var <- "pct_change_all"

### CORRESPONDING OPTIONS: variable types, labels, color schemes, and whether to reverse the color schemes

## VARIABLE UNITS/TYPES 
# percent, comma or dollar currently implemented -- if you have a different unit, ask for help or check out the scales package
map_var_type <- "percent"

## VARIABLE LABELS
# tip: "</br>" is HTML that creates a line break
map_var_label <- "Change in Low-Wage</br>Employment (Earnin)"

## COLOR SCHEMES
# "atlas" or "pink-purple" like atlas covariates
color_scheme_name <- "atlas"
raw_color_scheme <- c("#890024","#a94138","#c86e4f","#e19e6f","#f2cd97","#ffffc2","#ccdcb5","#98baa8","#6a9799","#447586","#195473")

# reverse the color scheme for each variable
reverse_color_scheme <- FALSE

# how many colors/deciles do you want?
number_of_colors = 10

# what opacity you want the shapes to be
opacity = 0.7

# number of extra pixels per pixel to have
resolution <- 4

# mapshot user agent
mapshot_useragent <- 
  paste("Mozilla/5.0",
        "(compatible; MSIE 10.6; Windows NT 6.1; Trident/5.0; InfoPath.2; ",
        "SLCC1; .NET CLR 3.0.4506.2152; .NET CLR 3.5.30729; .NET CLR 2.0.50727)",
        "3gpp-gba UNTRUSTED/1.0",
        sep = "")


#***********************************************************************************************************************************
## Loop through cities
#***********************************************************************************************************************************

# in case your zcta variable is named ZCTA
names(input_data) <- tolower(names(input_data))

for (city in cities){
  mapdata <- input_data %>% dplyr::filter(cityname == city)
  cityfn <- tolower(str_replace_all(city, " ", ""))
  
  
  #***********************************************************************************************************************************
  # Pull in the spatial data
  #***********************************************************************************************************************************
  
  zip_starts_with <- unique(str_extract(str_pad(mapdata$zcta, 5, "left", "0"), "^.{3}"))
  
  # loading in the census tract shapefile data
  shpdata <- zctas(starts_with = zip_starts_with, cb = TRUE, year = 2010)
  
  #***********************************************************************************************************************************
  # Join spatial data with your data
  #***********************************************************************************************************************************
  shpdata <- shpdata %>%
    mutate(zcta = as.numeric(as.character(ZCTA5)))
  
  shpdata <- geo_join(shpdata, mapdata, by = "zcta", how = "left")
  
  # check that we haven't dropped any zctas from your dataset
  if (length(mapdata$zcta) != length(intersect(mapdata$zcta, shpdata$zcta))){
    stop("ZCTAs have been dropped")
  }
  
  #***********************************************************************************************************************************
  # Set city-specific options
  #***********************************************************************************************************************************
  
  # set lat long and legend position based on city
  
  if(cityfn == "newyorkcity"){
    latlong <- c(-73.959128, 40.711723)
    legendpos <- "topleft"
    cityname <- "New York"
    output_name <- "NYC"
  } 
  if (cityfn == "sanfrancisco") {
    latlong <- c(-122.444719, 37.764026)
    legendpos <- "bottomleft"
    cityname <- "San Francisco"
    output_name <- "SF"
  } 
  if(cityfn == "chicago"){
    latlong <- c(-87.661129, 41.874609)
    legendpos <- "topright"
    cityname <- "Chicago"
    output_name <- "Chicago"
  } 
  
  #***********************************************************************************************************************************
  # Create color palette and legend labels
  #***********************************************************************************************************************************
  
  color_pal2 <- colorRampPalette(raw_color_scheme)
  
  quantiles <- unique(quantile(mapdata[[map_var]], 
                               prob = (1:(number_of_colors - 1)) / number_of_colors, na.rm = TRUE))
  
  color_palette <- color_pal2(length(quantiles))
  
  
  calc_quantiles <- quantiles / 100
  
  # number of decimal places in legend
  precision = 0.1
  
  
  formatter <- get(map_var_type)
  legend_labels <- paste0(formatter(calc_quantiles[-length(calc_quantiles)], accuracy = precision), 
                          " to ", 
                          formatter(calc_quantiles[-1], accuracy = precision))
  
  legend_labels <- c(paste0("< ", 
                            formatter(calc_quantiles[1], accuracy=precision)),
                     legend_labels,
                     paste0("> ", formatter(calc_quantiles[length(calc_quantiles)], accuracy = precision))) 
  
  quantile_palette <- colorBin(color_palette, mapdata[[map_var]], 
                               bins = c(-Inf, quantiles, Inf), 
                               reverse = reverse_color_scheme, 
                               right = FALSE)
  
  #***********************************************************************************************************************************
  # Produce the map
  #***********************************************************************************************************************************
  
  basemap <- 
    paste("https://api.mapbox.com/styles/v1/kaimatheson/",
          "ck21z1gw10ojz1cp9nohuo064/tiles/256/{z}/{x}/{y}@2x?",
          "access_token=pk.eyJ1Ijoia2FpbWF0aGVzb24iLCJhIjoiY2p6bXVvMGY3MGgwejNia3p5NGY4ZHNxbyJ9.4W1gdX9t8hveowSiVstyhQ",
          sep = "")
  
  basemap_labels <- 
    paste("https://api.mapbox.com/styles/v1/kaimatheson/",
          "ck21z38gi4nap1cmdw87c2h8n/tiles/256/{z}/{x}/{y}@2x?",
          "access_token=pk.eyJ1Ijoia2FpbWF0aGVzb24iLCJhIjoiY2p6bXVvMGY3MGgwejNia3p5NGY4ZHNxbyJ9.4W1gdX9t8hveowSiVstyhQ",
          sep = "")
  
  bm_attr <- "? <a href='https://www.mapbox.com/map-feedback/'>Mapbox</a>"
  
  mymap_og <- 
    leaflet(shpdata, 
            options = leafletOptions(zoomSnap = 0.01, zoomDelta = 0.01, attributionControl = FALSE)) %>% 
    addMapPane("background", zIndex = 410) %>%
    addMapPane("polygons", zIndex = 420) %>%
    addMapPane("labels", zIndex = 430) %>%
    # set basemap
    addTiles(urlTemplate = basemap, attribution = bm_attr,
             options = pathOptions(pane = "background")) %>% 
    addTiles(urlTemplate = basemap_labels, 
             options = pathOptions(pane="labels")) %>%
    addPolygons(color = ~quantile_palette(get(map_var)),
                weight = 1,
                fillOpacity = opacity,
                options = pathOptions(pane = "polygons")
    )
  
  for (image_size in c("papersize")){
    
    # set image width and zoom based on city and image size
    if (image_size == "normalsize"){
      vwidth = 1200
      
      if(cityfn == "newyorkcity"){
        zoom <- 11
      } 
      if (cityfn == "sanfrancisco") {
        zoom <- 10.45
      } 
      if(cityfn == "chicago"){
        zoom <- 12
      } 
    }
    if (image_size == "papersize"){
      vwidth = 700
      
      if(cityfn == "newyorkcity"){
        zoom <- 10.25
      } 
      if (cityfn == "sanfrancisco") {
        zoom <- 9.7
      } 
      if(cityfn == "chicago"){
        zoom <- 11.3
      } 
    }
    if (image_size == "nontechsize"){
      vwidth = 1500
      
      if(cityfn=="newyorkcity"){
        zoom <- 11
      } 
      if (cityfn == "sanfrancisco") {
        zoom <- 9.7
      } 
      if(cityfn=="chicago"){
        zoom <- 11.3
      }
    }
    
    mymap <- mymap_og %>% 
      setView(lat = latlong[2], lng = latlong[1], zoom = zoom) 
    
    mymap_legend <- mymap %>%
      addLegend(legendpos, pal = quantile_palette, values = ~get(map_var),
                title = map_var_label,
                opacity = 1, 
                labFormat = function(type, cuts) {
                  cuts = paste0(legend_labels)
                }, na.label = "No Data"
      ) 
    
    # save map as an image file
    vheight = 744/992 * vwidth
    
    filename <- paste("map_employment_byZIP_", output_name, file_ext, sep = "")
      
    mymap_legend %>% 
      mapshot(file = file.path(output_folder, filename),
              vwidth = vwidth, 
              vheight = vheight, 
              zoom = resolution,
              useragent = mapshot_useragent)
  }
}
